Figure 5

Heatmap of significant genes per cluster

hdWGCNA analysis

Load data and libraries

##################
# LOAD LIBRARIES #
##################
library(tidyverse)
library(Seurat)
library(SeuratObject)
library(tidyseurat)
library(hdWGCNA)
library(enrichR)
library(png)
library(cowplot)
library(patchwork)
library(openxlsx)

source("../../bin/spatial_visualization.R")
source("../../bin/plotting_functions.R")

#########
# PATHS #
#########
input_dir <- "../../results/06_DGE_condition_st_data/"
result_dir <- "./Figures/05/"
if( isFALSE(dir.exists(result_dir)) ) { dir.create(result_dir,recursive = TRUE) }
epi_clus <- "^5$|^6$|^7|^8" # res 0.7

ord <- c("Superficial", "Upper IM", "Lower IM", "Basal","1","4","0","3","2","9","10","11","12")
ord1 <- c("5", "6", "7", "8","1","4","0","3","2","9","10","11","12")
sample_id <- c("P020", "P045", "P050", "P057",
               "P008", "P031", "P044","P080", "P026", "P105", 
               "P001", "P004", "P014", "P018", "P087", "P118",
               "P021", "P024", "P067", "P081", "P117" ) %>% set_names()

#############
# LOAD DATA #
#############
# hdWGCNA
DATA <- readRDS(paste0("../../results/09_hdWGCNA/","hdWGCNA_3771DEGs_Seurat.RDS"))
modules <- read_csv(paste0("../../results/09_hdWGCNA/", "wgcna_3771DEGs_modules.csv"))

# DATA <- readRDS(paste0("../../results/09_hdWGCNA/",,"all_Clus_4000DEGs/","hdWGCNA_Seurat.RDS"))
# modules <- read_csv(paste0("../../results/09_hdWGCNA/all_Clus_4000DEGs/", "wgcna_all_Clus_modules.csv"))

# plot the dendrogram
png(paste0("./Figures/05/", "Spatial_hdWGCNA_dendrogram.png"), 
    width = 10, height = 3, units = "in", res = 300)
PlotDendrogram(DATA, main=NA, marAll = c(1, 4, 1, 0))
dev.off()

## quartz_off_screen 
##                 2

dend <- readPNG("./Figures/05/Spatial_hdWGCNA_dendrogram.png")
g <- grid::rasterGrob(dend, interpolate=TRUE) 

A <- ggplot() +
  annotation_custom(g, xmin=-Inf, xmax=Inf, ymin=-Inf, ymax=Inf) +
  theme_nothing() +
  theme(rect = element_blank(), # removes the box around the plot)
        plot.margin = unit(c(-0,0,0,0), "lines")) #t,r,b,l  

# dev.new(width=3.5, height=3, noRStudioGD = TRUE)
A

# get module eigengenes and gene-module assignment tables
MEs <-  DATA@misc[["vis"]][["MEs"]]
# add the MEs to the seurat metadata so we can plot it with Seurat functions
DATA@meta.data <- cbind(DATA@meta.data, MEs)

# plot with Seurat's DotPlot function
mods <- unique(modules$module)[unique(modules$module) != 'grey']
p <- c('groups', 'layers') %>%
  map(., ~DotPlot(DATA, features=rev(mods), group.by = .x, dot.min=0.1) )

# flip the x/y axes, rotate the axis labels, and change color scheme:
p <- map(p,~.x +
  coord_flip() +
  RotatedAxis() +
  scale_color_gradientn(colours = c('blue','grey95', 'red'), name = "Avg. Expression",
                       limits = c(-1.5, 2.5), oob = scales::squish,
                       values = scales::rescale(c(-1.5, 0, 2.5))) + 
    theme(axis.text = element_text(size = 10),
          legend.text = element_text(size = 9), legend.title =element_text(size = 10),
          axis.line = element_line(size = .4), # , colour = "#bebebe"
          axis.ticks = element_line(size = .4), # , colour = "#bebebe"
          legend.direction = "vertical", legend.box = "horizontal",
          legend.margin=margin(1,-1,1,6),
          plot.margin = unit(c(.7, -0, -1.2, -1), "lines")) + # t,r,b,l
    xlab('') + ylab('') + guides(col = guide_colourbar(barwidth = .3, barheight = 4 )) 
)

# combine plots
# dev.new(width=3.9, height=3.3, noRStudioGD = TRUE)
leg <- get_legend(p[[1]])
p1 <- plot_grid(p[[1]]+theme(legend.position = "none"),NULL, rel_widths = c(1, 1.3))
(B <- plot_grid(p1, p[[2]]+theme(legend.position = "none"),rel_heights = c(.75,1,1), ncol = 1) )

ggsave(filename=paste0("./Figures/05/", "dot-plot.png"),B,  width = 3.9, height = 3.3, bg = "white")

# relative expression level of each module
plot_list <- ModuleRadarPlot(
  arrange(DATA, layers),
  group.by = 'groups', combine = F, # ncol = 4, 
  #barcodes = seurat_obj@meta.data %>% subset(cell_type == 'INH') %>% rownames(),
  axis.label.size=3, group.line.width =  0.3, grid.line.width = 0.3, gridline.max.linetype = "dashed",
  grid.label.size=3
) 

# dev.new(width=5.5, height=2, noRStudioGD = TRUE)
(C <- wrap_plots(plot_list, ncol=4) & theme(title = element_text(size=8), plot.margin = unit(c(.1, -0, -1, -0), "lines")) )

# ggsave(filename=paste0("./Figures/05/", "Group_contribution.png"),C,  width = 10, height = 3, bg = "white")

plot_filt.fun <- function(DATA, gr = "L1"){
  DAT <- filter(DATA, groups == gr)
  DAT@misc[["vis"]][["MEs"]] <- DATA@misc[["vis"]][["MEs"]][colnames(DAT),]
  plot_list <- ModuleFeaturePlot(DAT, reduction = "umapharmony", features = "MEs", title =F) 
  return(plot_list)
}

mod <-  c("SM1", "SM2", "SM3", "SM4") # , "SM5", "SM6"
mod <- map(mod, ~plot_genes.fun(DATA, 
               gene = .x, 
               scale = F,
               mins = -20, maxs = 20,
               diverging = T,
               col = rev(c("#c41625","#dc4e43","#fa9975","#FDDBC7","#F7F7F7","#D1E5F0","#92C5DE","#4393C3","#2166AC")),
               point_size = .5,
               red="umapharmony", 
               lable = TRUE) + ggtitle(" ") + theme(plot.title = element_text(hjust = 0))) # size = 3

# dev.new(width=4, height=4, noRStudioGD = TRUE)
(D_1 <- map(seq_along(mod), ~mod[[.x]] + facet_wrap(~groups, ncol = 4)) %>% wrap_plots(., ncol=1) )

# ggsave(filename=paste0("./Figures/05/", "Modules_across_groups_UMAP.png"),D_1,  width = 10, height = 10, bg = "white")

# enrichr databases to test
dbs <- c('GO_Biological_Process_2021','KEGG_2021_Human','Transcription_Factor_PPIs')

# perform enrichment tests
DATA <- RunEnrichr( #map(dbs, ~
  DATA,
  dbs=dbs, # character vector of enrichr databases to test
  max_genes = Inf # number of genes per module to test. use max_genes = Inf to choose all genes!
)

## Uploading data to Enrichr... Done.
##   Querying GO_Biological_Process_2021... Done.
##   Querying KEGG_2021_Human... Done.
##   Querying Transcription_Factor_PPIs... Done.
## Parsing results... Done.
## Uploading data to Enrichr... Done.
##   Querying GO_Biological_Process_2021... Done.
##   Querying KEGG_2021_Human... Done.
##   Querying Transcription_Factor_PPIs... Done.
## Parsing results... Done.
## Uploading data to Enrichr... Done.
##   Querying GO_Biological_Process_2021... Done.
##   Querying KEGG_2021_Human... Done.
##   Querying Transcription_Factor_PPIs... Done.
## Parsing results... Done.
## Uploading data to Enrichr... Done.
##   Querying GO_Biological_Process_2021... Done.
##   Querying KEGG_2021_Human... Done.
##   Querying Transcription_Factor_PPIs... Done.
## Parsing results... Done.

# retrieve the output table
enrich_df <- GetEnrichrTable(DATA) %>%
  filter(Adjusted.P.value < 0.05) %>% 
  split(~db)

enrich_df %>%
  map(., ~split(.x, ~module)) %>%
  imap(., ~write.xlsx(.x, paste0(result_dir,"New_3771DEGs/", "Enrichment_",.y,".xlsx")) )

## $GO_Biological_Process_2021
## A Workbook object.
##  
## Worksheets:
##  Sheet 1: "SM1"
##  
## 
##  Sheet 2: "SM2"
##  
## 
##  Sheet 3: "SM3"
##  
## 
##  Sheet 4: "SM4"
##  
## 
##  
##  Worksheet write order: 1, 2, 3, 4
##  Active Sheet 1: "SM1" 
##  Position: 1
## 
## 
## $KEGG_2021_Human
## A Workbook object.
##  
## Worksheets:
##  Sheet 1: "SM1"
##  
## 
##  Sheet 2: "SM2"
##  
## 
##  Sheet 3: "SM3"
##  
## 
##  
##  Worksheet write order: 1, 2, 3
##  Active Sheet 1: "SM1" 
##  Position: 1
## 
## 
## $Transcription_Factor_PPIs
## A Workbook object.
##  
## Worksheets:
##  Sheet 1: "SM1"
##  
## 
##  Sheet 2: "SM2"
##  
## 
##  Sheet 3: "SM3"
##  
## 
##  Sheet 4: "SM4"
##  
## 
##  
##  Worksheet write order: 1, 2, 3, 4
##  Active Sheet 1: "SM1" 
##  Position: 1

# saveRDS(enrich_df, paste0("../../results/09_hdWGCNA/New_3771DEGs/", "Enrichment.RDS"))
# enrich_df <- readRDS(paste0("../../results/09_hdWGCNA/New_3771DEGs/", "Enrichment.RDS"))

#######################
# ENRICHMENT BARPLOT #
#######################
overlap.fun <- function(string){
  l <- str_split(string, pattern ="/")
  l <- map_dbl(l, ~as.numeric(.x[1])/as.numeric(.x[2]) )
  return(l)}

# str_match("cytoplasmic translation (GO:0002181)", "^(.+?)\\s\\((.+)\\)$")[2]
GeneRatio_plot.fun <- function(enrich_df, txt_size = 15, nr_path=3,
                               col=c("#ed968c","#f9d14a","#88a0dc","#e78429")){
  dot_df <- enrich_df %>%  # dot_df <- enrich_df$GO_Biological_Process_2021 %>%
    {if(grepl("GO",.$Term[[1]])){mutate(., "Term" = str_match(.$Term, "^(.+?)\\s\\((.+)\\)$")[,2], 
           "GOid" = str_match(.$Term, "^(.+?)\\s\\((.+)\\)$")[,3] )}else .} %>%
    filter(Adjusted.P.value < 0.05) %>%
    mutate("-log10(P-value)" = -log10(P.value)) %>%
    mutate(GeneRatio = overlap.fun(.$Overlap)*100) 
  
  dot_df <- dot_df %>%
    group_by(module) %>% 
    top_n(., -nr_path, Adjusted.P.value) 
  
  p <- ggplot(dot_df, aes(x = `-log10(P-value)`, y = fct_reorder(Term, `-log10(P-value)`), fill = module, col = module)) +
    geom_col(width = .05, show.legend = F) +
    geom_point(aes(size = GeneRatio)) + theme_classic() + 
    scale_fill_manual(values = col, aesthetics = c("fill", "colour")) +
    
    facet_wrap(~module, scales = 'free', ncol = 1) +
    #scale_x_continuous(expand = c(0, 1.5)) +
    coord_cartesian(clip = F) +
    scale_x_continuous(limits = function(x){c(0, +max(0.1, x))}) +
    scale_size_continuous(breaks = c(10, 50, 80)) +
    theme(axis.title.y = element_blank(),
          strip.text.x = element_text(hjust = 0.1, margin=margin(l=0)),
          strip.background = element_blank(),
          panel.spacing = unit(1, "lines"),
          axis.text.y = element_text(size = 10),
          panel.border = element_blank())
  return(p)
}

df <- bind_rows(enrich_df[[2]], enrich_df[[3]]) %>%
  filter(db == "Transcription_Factor_PPIs" & module == "SM4" | db == "KEGG_2021_Human")

# dev.new(width=3.5, height=5, noRStudioGD = TRUE)
(D_2 <- GeneRatio_plot.fun(df) + theme(legend.position = "none") )

# ggsave(paste0("./Figures/05/", "Enrichment.png"),D_2,  width = 3.8, height = 5, bg = "white")

######################
# COMBINE ALL PANELS #
######################
# dev.new(width=7, height=8, noRStudioGD = TRUE) 
A_C <- plot_grid(A, C, ncol = 1, rel_widths = c(1, .4))
A_B_C <- plot_grid(A_C, B, ncol = 2)
D <- plot_grid(D_1, D_2)
(Figure5 <- plot_grid(A_B_C, D, ncol = 1, rel_heights = c(0.33, .66)) )